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Abstract 



The Turing instability paradigm is revisited in the context of a multispecies 
diffusion scheme derived from a self-consistent microscopic formulation. The 
" ^ . analysis is developed with reference to the case of two species. These latter 

share the same spatial reservoir and experience a degree of mutual interfer- 
p I ence due to the competition for the available resources. Turing instability can 

set in for all ratios of the main diffusivities, also when the (isolated) activa- 
tor diffuses faster then the (isolated) inhibitor. This conclusion, at odd with 
the conventional vision, is here exemplified for the Brusselator model and 
ultimately stems from having assumed a generalized model of multispecies 
diffusion, fully anchored to first principles, which also holds under crowded 
conditions. 

Keywords: Turing instability. Stochastic processes, Reaction-diffusion 
systems. Cross-diffusion systems 



1. Introduction 

Turing instability is one of the reference mechanisms for pattern forma- 
tion in nature [H, [ij . The Turing idea applies to a large gallery of phenomena 
that can be modelled via reaction-diffusion equations These latter are 
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mathematical models that describe the dynamical evolution of distinct fam- 
ihes of constituents, mutually coupled and freely diffusing in the embedding 
medium. Diffusion can seed the instability by perturbing the mean-field ho- 
mogeneous state, through an activator-inhibitor mechanism, and so yielding 
the emergence of patched, non homogeneous in space, density distributions. 
The most intriguing applications of the Turing paradigm are encountered in 
the context of morphogenesis, the branch of embryology which studies the 
development of patterns and forms in biology. The realm of application of 
the Turirig ideas encompasses however different fields, ranging from chem- 
istry [19|, l20|, l2l| to biology Il8| , passing through physics [10| , 
where large communities of homologous elements evolve and interact. 

According to the classical viewpoint, however, the diffusion coefficient of 
the inhibitor species has to be larger than that of the activator, for the pat- 
terns to eventually develop. This is a strict mathematical constraint which 
is not always met in e.g. contexts of biological relevance 23|, |22|], and which 



limits the possibility of establishing a quantitative match between theory and 
empirical data. Spatially extended systems made of interacting species shar- 
ing similar diffusivities can indeed display self-organized patched patterns, 
an observation that still calls for a sound interpretative scenario, beyond the 
classical Turing mechanisms 

One viable strategy to possibly reconcile theory and observations has 



been explored in [ll|] and [12|. In these studies, the authors considered the 
spontaneous emergence of persistent spatial patterns as mediated by the de- 
mographic endogenous noise, stemming from the intimate discreteness of the 
scrutinized system. The intrinsic noise translates into a systematic enlarge- 
ment of the parameter region yielding the Turing order, when compared to 
the corresponding domain predicted within the deterministic linear stability 
analysis. It is however unclear at present whether experimentally recorded 
patterns bear the imprint of the stochasticity, a possibility that deserves to 
be further challenged in the future. 

Alternatively, and to bridge the gap with the experiments, the Turing 
instability concept has been applied to generalized reaction-diffusion equa- 
tions. These latter account for cross diffusion terms which are hypothesized 
to exist on purely heuristic grounds or by invoking the phenomenological 



theory of linear non-equilibrium thermodynamics [13|, |16|, ll7[. Diagonal and 



off-diagonal coefficients of the diffusion matrix are not linked to any micro- 
scopic representation of the examined dynamics and are hence treated as 
free parameters of the model. In IJ] the authors quantify the impact of 
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cross terms on the Turing bifurcation, showing e.g that spatial order can 
materiahze also if the inhibitor's diffusion ability is less pronounced than the 
activator's one. 

Starting from this setting, the aims of this paper are twofold. On the one 
side, we shall elaborate on a microscopic theory of multispecies diffusion, fully 
justified from first principles. The theory here derived is specifically targeted 



to the two species case study and extends beyond the formulation of [15 
On the other side, and with reference to the Brusselator model, we will show 
that Turing patterns can take place for any ratio of the main diffusivities. 
In doing so we will cast the conclusions of 14| into a descriptive framework 
of broad applied and fundamental interest, where the key cross diffusion 
ingredients are not simply guessed a priori but rigorously obtained via a 
self-consistent derivation anchored to the microscopic world. Working in the 
context of a reference case study, the Brusselator model, we shall also perform 
numerical simulations based on both the underlying stochastic picture and 
the idealized mean-field formulation to elaborate on the robustness of the 
observed patterns. 

In the following we briefly discuss the derivation of the model, focusing 
on the specific case where two species are supposed to diffuse, sharing the 
same spatial reservoir. 



2. The model 

Consider a generic microscopic system bound to occupy a given volume 
of a (i— dimensional space. Assume the volume to be partitioned into a large 
number Q of small hypercubic patches, each of linear size /. Each mesoscopic 
cell, labelled by i, is characterized by a finite carrying capacity: it can host 
up to N particles, namely nf of type A, nf of type B, and Vi = N — nf — nf 
vacancies, hereafter denoted by V . In general, the species will also interact, 
as dictated by specific reaction terms. Let us start by solely focusing on 
the diffusion part, silencing any direct interaction among elementary con- 
stituents. As we shall remark, there exists an indirect degree of coupling 
that results from the competition for the available spatial resources. In prac- 
tice, the mobility of the particles is balked if the neighbouring patches have 
no vacancies. Particles may jump into a nearest-neighbour patch, only if 
there is a vacancy to be eventually filled. This mechanism translates into the 
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following chemical equation 



(1) 



where i and j label nearest-neighbour patches. Here, Ai and Bi identify the 
particles that belong to cell i. Vi labels instead the empties that are hosted in 
patch i. The parameters /i"^ and fi^ stand for the associated reaction rates. 
Similar reactions control the migration from cell j towards cell i. 

In addition, and extending beyond the scheme proposed in 15j, we imag- 
ine the following reactions to hold: 



A, + Bj 
Ai + Bi 



A, + Bi 
Ai + B, 



(2) 



which in practice account for the possibility that elements Ai (resp. Aj) and 
Bj (resp. Bi) swap their actual positions. 

The state of the system is then specified by the number of A and B 
particles in each patch, the number of vacancies following from a straightfor- 
ward normalization condition. Introduce the vector n = (rii, . . . , n^), where 



Hi 



71-, n. 



). The quantity T(n'|n) represents the rate of transition from 
state n, to another state n', compatible with the former. The transition rates 
associated with the migration between nearest-neighbour, see Eqs. ([1]), take 
the form 



Tfn^") - 1. 



+ l\nl \ rV- \ 



zVt N 



N 



A,B, 



(3) 



where we have made explicit in the components that are affected by 

the reactions. As discussed in isl 



the factor N — — , 



reflects the 

natural request of a finite capacity, and will eventually yield a macroscopic 
modification of the Pick's law of diffusion. Moreover, chemical equations (E]) 
result in the following transition rates: 



T{nf - 1, nf + 1, nf + 1, nf - , nf, nf, nf) 
Tint + 1, nf - 1, nf - 1, nf + l\nf, nf, nf , nf) 



a nf rij 



a nfnf 
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The process here imagined is Markov, and the probabihty P(n, t) to ob- 
serve the system in state n at time t is ruled by the master equation 

= J2 [^(nln')P(n', t) - T(n'|n)P(n, t)] , (5) 

where the allowed transitions depend on the state of the system via the 
above relations. Starting from this microscopic, hence inherently stochastic 
picture, one can derive a self-consistent deterministic formulation, which 
exactly holds in the continuum limit. Mathematically, one needs to obtain the 
dynamical equations that govern the time evolution of the ensemble averages 
(nf-) and (nf ). To this end, multiply first the master Eq. by n", with a = 
A, B, and sum over all n. After an algebraic manipulation which necessitates 
shifting some of the sums by ±1, one eventually gets 



dt 



= J2[{T{nt^ + l,nf-l\nt\nf)) (6) 

+ {T{nt + l,nf-l, nf - 1, nf + l|nf , nf, nf , nf)) 
-{Tin^^-l,nf + l\n'f'\nf)) 
- {Tint - 1, n/ + 1, nf + 1, nf - 1 |nf , nf, nf, nf)) 

where the notation Cleans that we are summing over all patches j which 
are nearest-neighbours of patch i. The averages in Eq. ([6]) are performed 
explicitly by recalling the expression for the transition rates as given in Eqs. 
([3]) and (jl]). Replace then the averages of products by the products of av- 
erages, an operation that proves exact in the continuum limit N ^ oo. By 



introducing the continuum concentration {(pA,B)i = lini7v_j.oo , rescaling 
time by a factor of NQ and taking the size of the patches to zero one finally 
get^ 



dt 
dt 



^Use has been made of the discrete Laplacian operator Afi — [2/ z)J2j!zi{fj ^ li)^ 
which then turns into the continuum operator V when sending to zero the size of the 
patch and scahng the rates fji^'^ and a appropriately. 
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wher^ -Dii^22 1'^I^a,b and -Di2,2i — ^ P{i^a,b — «)• The above system of 
partial differential equations for the concentration 0^ and 0^ is a slightly 



modified version of the one derived in [15|, this latter being formally recovered 



when setting a to zero. In the generalized context here considered, the cross 
diffusion coefficients D12 and D21 are different, specifically smaller, than the 
corresponding mean diffusivities Dn and D22- We emphasize again that the 
crossed, nonlinear contributions ±(0^^bV^0b,a — 4'b,a^'^4'a,b) stem directly 
from the imposed finite carrying capacity and, as such, have a specific, fully 
justified, microscopic origin. The diffusive fluxes that drive the changes in 
the concentrations (pA and (ps can be written as: 

J^s = -^210i?V0A - ^22 (^1 - V0B (8) 

It is interesting to notice that relations ([8]) enable us to make contact 
with the field of linear non-equilibrium thermodynamics (LNET), a branch 
of statistical physics which defines the general framework for the macroscopic 
description of e.g. transport processes. One of the central features of LNET 
is the relation between the forces, which cause the state of the system to 



change, and the fluxes, which are the result of these changes [13|. Within 
the formalism of LNET the fluxes and J,^^ that rule the diffusion of the 
two species 0^ and 4'b are linearly related to the forces, the gradients of 
the respective concentrations. The quantities that establish the formal link 
between forces and fluxes are the celebrated Onsager coefficients, postulated 
on pure heuristic grounds. Interestingly, Eqs. (|8]) provide a self-consistent 
derivation for the Onsager coefficients, that enters the generalized Pick's 
scenario here depicted. 

Define $ = (0a5 0b) and J = (J</,^, J,Pb)- Then Eqs. ([7]) can be written 
in the compact form: 

_ = -VJ = VD($)V$ (9) 

"^From the above expressions, one derives the consistency conditions > <^ and fis > 

rv 
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where the 2x2 matrix D reads: 



D($) 



Du (l - 



) 



) 



) 



A stringent constraint from thermodynamics is that all eigenvalues of 
the diffusion matrix D are real and positive. This in turn corresponds to 
requiring tr(D) > and det(D) > 0. A straightforward calculation yields: 



where AD = Du — D12 = -D22 — ^21- By definition AD > 0. Moreover, 
(pA and (pB are both positive and smaller than one. Hence, tr(D) > and 
det(D) > 0, a result that points to the consistency of the proposed formula- 
tion. 

3. The region of Turing order 

Having derived a plausible macroscopic description for the two compo- 
nents diffusion process, we can now move on by allowing the involved species 
to interact and consequently consider in the mathematical model the cor- 
responding reaction terms. As an important remark, we notice that these 
latter can be also obtained as follows the above, rather general, approach 
that bridges micro and macro realms. First, one need to resolve the interac- 
tions among individual constituents, by translating into chemical equations 
the microscopic processes implicated. These include cooperation and com- 
petition effects, as well as the indirect interferences stemming from the finite 
carrying capacity that we have imposed in each mesoscopic patch. Then, 
one can recover the deterministic equations for the global concentrations, by 
operating in the continuum system size limit. In general, Eq. is modified 
into: 



where F = {fA{(pA,(pB),fB{4>A,(pB))- As we have anticipated, the interest of 
this generalized formulation, resides in that it allows for Turing like patterns 
in a region of the parameter space that is instead forbidden when conventional 
reaction-diffusion systems are considered. The novelty of the proposed for- 
mulation has to do with the presence of specific cross diffusion terms, which 



tr(D) 
det(D) 



Du{l - 4>b) + D22{1 - (Pa) + AD{(Pa + <Pb) 
^11^22(1 -<t>A- <Pb) + AD{Du(pA + D22<Pb) 




(10) 
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follow a sound physical request, and add to the classical Laplacians, signature 
of Fickean diffusion. 

Let (t)A,4'B be the steady state solution of the homogeneous (aspatial) 
system, namely fA{,4>A,4>B) = fB{4>A,(t>B) = 0. The fixed point is linearly 
stable if the Jacobian matrix A 

A = 



dfA dfA 

dJ>A dJ>B 

dfs dfB 



has positive determinant and negative trace. It is worth stressing that the 
derivatives in matrix A are evaluated at the homogeneous fixed point. Back 
to the complete model, a spatial perturbation superposed to the homogeneous 
fixed point can get unstable if specific conditions are met. Such conditions, 
inspired to the seminal work by Turing, are hereafter derived via a linear 
stability analysis. Define rj = ^ — ^ and proceed with a linearization of Eq. 
( ITOl) to eventually obtain: 

at 

Going to Fourier space one gets: 

§ = A-m (11) 

where A*(A;) = A($) — A;^D($). By characterizing the eigenvalues of the 
matrix A*, one can determine whether a perturbation to the homogeneous 
solution can yield patterns formation. In particular, if one of the eigenvalues 
admits a positive real part for some values of k, then a spatially modulated 
instability develops. The growth of the perturbation as seeded by the linear 
instability will saturate due to the non linearities and eventually results in a 
characteristic pattern associated to the unstable mode k. Steady patterns of 
the Turing type require in addition that the imaginary part of the eigenvalues 
associated to the unstable mode are zero. In formulae, the Turing instability 
sets in if there exists a k such that tr(A*(A;)) < and det(A*(A;)) < 0. 
These latter conditions are to be imposed, jointly with the request of a stable 
homogeneous fixed point (tr(A) < 0, det(A) > 0), to identify the parameters' 
values that drive the instability. Alternatively, one can obtain a set of explicit 



conditions following the procedure outlined below, and adapted from [12 
The eigenfunctions of the Laplacian operator are: 

(V2 + Wfc(r) = 0, 
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and we write the solution to Eq. fllip in the form: 



x(t,r) = ^e^* au Wfc(r) 



(12) 



By substituting the ansatz (fT^ into Eq. (ITT]) yields: 

The above system admits a solution if the matrix A — /c^ D — Al is singular, 
i.e.: 

det(A - A;^ D - Al) = 0. (13) 

The solutions A(A;) of ( TT^ can be interpreted as dispersion relations. If 
at least one of the two solutions displays a positive real part, the mode is 
unstable, and the dynamics drives the system towards a non-homogeneous 
configuration in response to the initial perturbation. Introduce the auxiliary 
quantity V defined as: 



dfB dfA 
d(f)B d(j)A 



dfA dfB 

d(f)A d(j)A 



D 



B 



12" 



D 



dfA 



21" 



^B 



^B 



(14) 



Then a straightforward calculation results in the following compact con- 
ditions for the instability to occur: 



r > 

> ADiiD 



(15) 



22 



D 



12 



D 



21 



D 



22 



B det(A) 



together with tr(A) < and det(A) > 0. 

For demonstrative purposes we now specialize on a particular case study 
and trace out in the parameters' plane, the domain that corresponds to the 
Turing instability. Our choice is to work with the Brusselator mode0 which 
implies setting = -{h + d)(t)A + a(l - 4>a - 4>b) + c(f)\(f)B and = b(f)A - 



^The term a(l — (j>A ~ (f's) reflects the presence of the finite carrying capacity, as 
discussed in [l^ . Similar conclusions hold however if the diluted limit is performed, just 
in the reaction terms, hence replacing a(l — (j)A — (t>B) with a. 
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c4>\4>B- Species A plays now the role of the activator, while B stands for the 
inhibitor. Results of the analysis are reported in left panels of Fig. [H where 
the region of interest is singled out in the plane {b, c), for different choices of 
AD. Turing patterns are predicted to occur for D22/-D11 < 1, at odd with 
what happens in the conventional scenario where standard Pick's diffusion 
is assumed to hold (see below). The right panels report the results of direct 
simulations and confirm the presence of macroscopically organized patterns 
in a region of the parameters space that is made classically inaccessible by the 
aforementioned, stringent condition D22 > -Dii The simulations refers to the 
choice D22/D11 = 0.7. These observations are general and similar conclusions 
can be drawn assuming other reactions schemes of the inhibitor/activator 
type, different from the Brusselator model. 

It is now instructive to elaborate on a simple interpretation of the above 
result. Let us start by briefly revisiting the necessary conditions for the 
classical Turing instability to occur, namely: 

tT{A) = dfA/d<f)A + dfB/d<j)B < (16) 
DndfB/d<l)B + D22dfA/d<j)A > 0. 

Both conditions can be simultaneously matched, only if the diagonal elements 
of the Jacobian matrix A have opposite signs. For the sake of clarity, let us 
assum^ that: 

d(t)A 0(l)B 

Hence, species A activates its own production, while species B has a self- 
inhibitory feedback. Requiring tr(A) < implies imposing |||^| > 
which, by making use of the second of ( |T6l) . readily translates into the nec- 
essary condition 

D22 \dfB/d(j)B\ , . 

Ai dfA/d<PA ^ ^ 

As already mentioned, the inhibitor must diffuse faster than the activator 
(when the two species are evolved in separate containers) for the conventional 
Turing pattern to occur: the system has to accommodate for two competing 



^This is indeed the case for the Brusselator model. For c sufficiently large, see also 
panels (a) and (c) of Fig. [U we have in fact ~ 2c4)a4'b > and = ~c$a < 0- 
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processes, a short-range activation and long-range inhibition. Starting from 
this setting we can adapt the above reasoning to the generahzed case study 
where cross diffusion terms are also present. To this end, and to keep the 
notation light, we shall solely consider the limiting case with AD = 0. Similar 
conclusions hold when AD ^ 0. The second of relations f lT6|) is now replaced 
by the condition F > (see Eq. f[T^ ). which can be cast in the form: 



D 



11 



df. 



B 



+ D 



22 



dcpB 

dfA 

d(f)A 



1 - 



B 



d<pA 
dfA 
'd(f)B 



(18) 



> 



when Dii = D12 and D22 = -D21. To proceed in the discussion we note that 
the elements that enter the square brackets have dimension of the inverse 
of time. Assume ^1 — 0^ j — 4>a^^ to be negative as it is reasonable 
to hypothesize if (i) the correction term that scales to the number densities 
0^ is sufficiently small, or conversely if (ii) we require dfs/dcpA > (i.e. 
the first species stimulates with a positive feedback the other). Under these 
conditions, one can then introduce the characteristic time scale tb associated 
to the reaction dynamics of species B, defined as: 



Tb 



df 



B 



90 



B 



B + 



df 



B 



d<P, 



(19) 



Similarly, for species A, we have: 

'dfA 



ta 



>B 



dh 



B 



-1 



(20) 



assuming dfA/ dip a to control the sign in the above expression, or alterna- 
tively imposing dfA/dcpB < (i.e. the second species acts with a negative 
feedback on the first one). The necessary condition fllSp for the generalized 
Turing instability to occur takes the form: 



'A 



taDu < tbD 



22 



where we have introduced two characteristic length scales, respectively I a, 
Ib, associated to the reactive dynamics of species A and B. In practice, also 
when D22 < Dii, spatially organized patterns can develop in the generalized 
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reaction diffusion scheme provided the activator has a shorter hfe time, than 
the inhibitor. In formulae, = T5-D22/-D11 < tb ■ In practical terms, the 
competition for the microscopic spatial resources modifies the time scales as- 
sociated to the reactions processes and induces a self-consistent long-range 
effect that enlarges the region of influence of the (isolated) inhibitors, also 
when the microscopic diffusion of the (isolated) activator is assumed to be 
faster. The crossed terms in the diffusion matrix determine a non trivial 
modification of the underlying characteristic times, which are now also sen- 
sitive to the off-diagonal elements of the Jacobian matrix. In the diluted 
limit in fact, ta — )■ t^^ = {dfA/d(f)A)~^ and tb — )■ r^*' = |c^/b/c?0_b|~^ and 
one is brought back to the standard, stringent condition (fT7|) . In Fig. |2]the 
ratio tb/ta is displayed for the Brusselator model, inside the Turing region, 
as a function of the chemical parameter b. Different curves refer to distinct 
choices of c, while the other parameters are set to the values of Fig. [1^, with 
AD = 0. As expected, t^/ta > 1 a condition that eventually yields the gen- 
eralized Turing patterns as described above. Conversely, and as pictured in 
the small inset, Tb^ /rj^^ < 1. Hence, since D22 < -Dn, Turing patterns cannot 
manifest via the classical pathway, which applies to diluted conditions. 

The remaining part of this section is devoted to discussing the robust- 
ness of the patterns depicted in Fig. [1] (panels (b) and (d)), and obtained 
upon integration of the governing system of partial differential equations. 
It should be emphasized however that the model of multispecies diffusion 
here considered is stochastic in nature. It is therefore interesting to further 
elaborate on the contributions played by finite size effects, associated to the 
graininess of the system, and hence deliberately neglected under the idealized 
deterministic representation of the dynamics. To this aim, one can carry out 
stochastic simulations, based on the Gillespie algorithm 2^, which produces 



realizations of the dynamics formally equivalent to those obtained from the 
governing master equation ([5]). We have here chosen to operate for the pa- 
rameter setting of Fig. [T]d and the results of our analysis are reported in 
Fig. 131 If the number of elements is sufficiently large (A^ = 3000, in the 
left panel of Fig. [3]) the patterns appear robust and resemble those recorded 
when operating in the framework of the deterministic picture. However, if 
the total number of microscopic individuals is reduced (A^ = 300, in the right 
panel of Fig. [3]) the patterns are less distinct and eventually fade away. De- 
mographic fluctuations ultimately destroy the self-organized spatial patterns, 
relic of Turing instability, and the system evolves towards an asymptotically 
stable homogeneous solution. The lifetime of the metastable non homoge- 
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neous patterns increases with the system size and formally diverges in the 
thermodynamic limit — > oo. Waiting for a sufficiently large time, also the 
apparently stable density structures as displayed in Fig. [3^1 are expected to 
coalesce and smear out. In other words, and intriguingly enough, the two 
limits for N oo and t — ?■ oo do not commute. If the system size limit is 
taken before the infinite time limit, the dynamics is permanently frozen into 
a stationary non homogeneous configuration, the spatially ordered Turing 
patterns. Conversely, the system is attracted towards a stable homogeneous 
equilibrium, due of the microscopic mixing that is seeded by the finite size 
fiuctuations. Clearly the time of homogeneization can be extremely long, 
when compared to the finite time window of the experimental observation. 
In this respect, the metastable spatially extended patterns are possibly the 
solely regimes to be accessible to direct measures. This observation shares 
many similarities with the phenomenon of Quasi-Stationary States, so far 



associated to the long range nature of the two-body interaction [26|, 125 
These findings, as well as the analysis of [2?!], can possibly shed new light 
onto the emergence of the Quasi-Stationary States, beyond the domain of 
applications for which they have been reported to occur. As a side remark, 
it is worth emphasising that similar conclusions hold when considering the 
diluted limit, i.e. when neglecting the role of a finite carrying capacity and 
the competition for the finite spatial resources that eventually yield the gen- 
eralized cross diffusion terms here considered. 



4. Conclusions 

Summing up, Turing patterns can develop for virtually any ratio of the 
main diffusivities in a multispecies setting. This striking effect originates 
from the generalized diffusion theory that is here assumed to hold and that 
builds on the scheme discussed in |15j]. Because of the competition for the 
available resources, a modified (deterministic) diffusive behaviour is recov- 
ered: cross diffusive terms appear which links multiple diffusing communities 
and which add to the standard Laplacian terms, relic of Pick's law. The fact 
that Turing like patterns are possible for e.g. equal diffusivities of the species 
involvec|§, as follows a sound dynamical mechanism, constitutes an intriguing 



^Notice that the authors of 12| failed to reahze that accounting for cross diffusion 
terms of the type derived in 15[ could result in an extension of the Turing mechanism to 
regions where D22 ^ Dn. 
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observation that hold promise to eventually reconcile theory and experimen- 
tal evidences. The investigated setting applies in particular to multispecies 
systems that evolve in a crowded environment, as it happens for instance 
inside the cells where different families of proteins and other biomolecular 
actors are populating a densely packed medium. It is interesting to notice 
that the stochastic fluctuations, endogenous to the scrutinized system in its 
discrete version, eventually destroy the patterns, that are instead deemed to 
be stable according to the idealized deterministic viewpoint. The lifetime of 
the metastable patched patterns increases however with the size of the sys- 
tem, in striking analogy with what has been observed for the so called Quasi- 
Stationary States, out-of-equilibrium regimes observed in systems subject to 
long-range interactions. For large enough N, the homogeneization as seeded 
by fluctuations is progressively delayed and eventually prevented in the con- 
tinuum limit N ^ oo. 
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b Iterations of the deterministic simulation 

(c) (d) 

Figure 1: Panels (a) and (c): the boundaries of the region of Turing instabihty are traced 
in the plane (&, c), for D22/D11 — 1 (panel (a)) and D22/D11 — 0.7 (panel (c)). The 
calculated domains refer to the Brusselator model with non Fickean diffusion, as explained 
in the main text. The solid line, which encloses regions I and II, stands for AD = 0, while 
the dashed line delimits region I, where the condition AD = 0.1 applies. The other 
parameters are set as a = 5, 0? = 3. Panels (b) and (d): the time evolution of the 
concentration (jjA, as revealed by direct numerical simulations. In both small 
perturbation is superposed at < = to the (non trivial) stable homogeneous fixed point of 
the Brussellator, namely 4>a = {a. + \/ — ^ah(a + d)/c)/2/{a + d), (j)B = hlc/4>A- Here, 
Dii = 1.0, D22 = 0.7, h = 21.71, c = 139, a = 5, d = 3. The upper right figure, panel (b), 
refers to AD = 0, the lower right, panel (d), to AD = 0.1. In the simulations we have 
assumed a symmetric box [— L,L], with L = 10. The box is discretized in 200, uniformly 
spaced, mesh points. The simulations are run by employing an explicit Euler scheme with 
time step equal to 0.0001. The density in each cell of the mesh is displayed in the vertical 
axis, while the horizontal axis refers to the number of iterations. 
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Figure 2: Main figure: the ratio tb/ta is plotted for the Brusselator model, inside the 
Turing region, as a function of the chemical parameter 6, for different choices of c. From 
left to right, c = 139, 139.4, 139.8. The other parameters are set as in Fig. [1^, with 
/S.D = 0. tb and ta follow respectively Eqs. and (^0)) and quantify the time scales of 
the reactive processes, within the framework of the generalized reaction diffusion scheme. 
As expected, the existence of a region of Turing order, as revealed in Fig [1^, implies 
Tb > TA- In the inset, the ratio of the time scales t^*'/t^*' obtained in the diluted limit is 
reported and proven to be smaller than unit. 




Figure 3: Time evolution of the discrete concentration ua/N, as it results from a direct 
integration of the stochastic Brusselator model. The simulations follows the Gillespie 



algorithm [2J|. Parameters refer to region II of Fig. [T}), namely Du — 1.0, D22 — 0.7, 



AZ? = 0, a = 5, d = 3, 6 = 21.71, c = 139. In panel (a): N = 3000, while in panel (b) 
N = 300. Demographic fluctuations destroy the deterministic patterns which are hence 
interpreted as a metastable regime of the finite N stochastic dynamics. 
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